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We present a method utihzing the continuity equation for the condensate density to make predic- 
tions of the precessional frequency of single off-axis vortices and of vortex arrays in Bose-Einstein 
condensates at finite temperature. We also present an orthogonalized Hartree-Fock-Bogoliubov 
(HFB) formalism. We solve the continuity equation for the condensate density self-consistently 
with the orthogonalized HFB equations, and find stationary solutions in the frame rotating at this 
frequency. As an example of the utility of this formalism we obtain time-independent solutions for 
quasi-two-dimensional rotating systems in the co-rotating frame. We compare these results with 
time-dependent predictions where we simulate stirring of the condensate. 
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I. INTRODUCTION 

One of the characteristics of a superfluid is the quan- 
tization of vortices that are found when sufficient angu- 
lar momentum is present in the system. Vortices, as a 
signature of superfluidity, are therefore of great theoret- 
ical and experimental interest, and the reader is referred 
to the recent review article by Fetter [ij. Our model 
is based on the Hartree-Fock-Bogoliubov (HFB) formal- 
ism We demonstrate how the continuity equation 
for the condensate density, solved self-consistently with 
a set of time-independent orthogonalized Hartree-Fock- 
Bogoliubov (HFB) equations in the frame rotating at 
the precessional frequency, can be used to make a pri- 
ori predictions of the precessional frequencies of vortices 
in Bose-Einstein condensates (BECs). By introducing a 
set of modified basis functions incorporating the vortex 
positions one is able to do this not only for single ofF-axis 
vortices, but also for the multiple vortex case. In order 
to perform these calculations correctly, it is necessary 
that the condensate and thermal populations be mutu- 
ally orthogonal. We present an orthogonal HFB formal- 
ism in which this condition holds, and for which we show 
the existence of a zero-energy eigenvalue in the time- 
independent case, in contrast with the standard HFB for- 
malism. As an illustration we model a two-dimensional 
BEG system, establishing the dependency of the preces- 
sional frequency on the lattice parameter a and on 
the temperature T for triangular and hexagonal vortex 
arrays. We also show the relationship of with T for 
the case of two vortices and for triangular and hexago- 
nal vortex arrays with three and seven vortices, respec- 
tively, having lattice parameter a = 3 harmonic oscilla- 
tor units, establishing the existence of an upper bound 
for the precessional frequency as a function of the num- 
ber of vortices. We obtain qualitative agreement with 
the areal density approximation d, 01 (see also Ref. Q) 
for the hexagonal vortex array. We verify the validity of 
these predictions in a series of finite t emp erature simu- 
lations using time-dependent HFB [§, M- We use one 
or more Gaussian stirrers to impart angular momentum 



to the BEG, and establish a critical stirring frequency 
required for the creation of vortices corresponding to a 
local stirring velocity just in excess of the local Landau 
critical velocity. We also verify that the axial component 
of the angular momentum is conserved when the trapping 
potential is axially symmetric, thus satisfying Noether's 
theorem. We show that angular momentum is lost when 
this symmetry is broken, leading to the decay of vortices. 



II. FORMALISM 

We consider a BEG system in the frame rotating with 
angular frequency CI with grand-canonical Bose Hamil- 
tonian given by 



where /ia(r) is the single-particle Hamiltonian 



(1) 



Mr) = -^V 



ihn- {v xV) + Vt{v), (2) 



and 



(3) 



with fls the s-wave scattering length. We wish to ob- 
tain an HFB formalism such that the condensate and 
thermal modes are orthogonal. We proceed by splitting 
the Bose field operator ^/'(r,t) into a coherent part rep- 
resented by the condensate field operator <l>(r,i), and an 
incoherent part represented by the fluctuation operator 
fi{r,t), writing i/'(r,t) = $(r, + fi{r,t) with $(r,t) = 
(j)(r,t)dc{t), with dc(t) the annihilation operator for the 
condensate, and (/)(r, t) a condensate wave function satis- 
fying the normalization condition / dr \(l>(r, t)\ = 1. The 
fluctuation operator fj{r,t) is defined as [ll!] 'fj{r,t) = 



2 



'0(r, i) — 4'{v,t) J drcj)* (r,t)il;{r,t). It can then be shown 
that the orthogonahty condition 



dr(l)*{r,t)ri{r,t) = 



(4) 



holds, so the condensate and thermal populations are 
orthogonal, as required. We also find from the definition 
of fj that dc{t) = J drcj)* {r,t)il;{r,t). In addition $ and i) 
satisfy the commutation relations 



<i>(r,i),Vt(i.',i) 



(r,t)0*(r',t). 



= 



(5) 



and 

r)(r, t), (r', t)] = 5{r - r') - 0(r, t)4>* (r', t) = Q{r, r', t), 



= 



(6) 



respectively. We use the Bogoliubov transformation 



to diagonalize the grand-canonical Hamiltonian (as with 
standard HFB), where we assume that the operators dk 
obey the Bosonic commutation relations 



dk,d] 



= Ski and [dk,di] = 



= 0. 



(8) 



We use the respective Heisenberg equations of motion for 
the operators $ and fj and the commutation relations 



dg,ip{r,t) =v*{r,t), 
a^,V't(r,t) =u*{r,t) 



(9) 



for the quasi-particle annihilation and creation operators 
dq and cij to derive the equations 



drcj)* 



9$ /- -\ - 

ih— + l^hn - M + -0 



(10) 



for the condensate operator $, and 



ih 
and 



duq{r,t) 



dt 



dv' \ l{y, r', i)i*9(r', t) + M (r, r', t)vq{Y' , t) 

(11) 



dt 



dv' [L(r, r', (r', + M(r, r', i)^i;(r', i) 

(12) ■ 

for the quasi-particle amplitudes Uq and Vg, where we 
have defined the operators 



L(r,r',i) = Q[v,v',t)[h{v')-^, + 2g^^v',t)^{v',t) 
M{v,v',t) = Q(r,r',t)(g^(r',t)V-(r',t)) 

(13) 



for notational convenience. We note that the quantity in 
square brackets in Eq. (fTO| is orthogonal to 0*. Hence 
we can re- write (ITTO in the form 



ih 



= (^o - + 5^^!-, t)^(r, t)) ^(r, t) - i(r, 

(14) 

where A(r, t) is a quantity orthogonal to 0, i.e. 
J (ir0*yi(r,t) = chosen such that angular and linear 
momentum conservation hold. Using the mean-field ap- 
proximations dc4> — > fj^fl — > (fff]) = n, fifj {fji)) = 

ih , fj^fj^ {v^V^) = "^*: and 'ip'^ 'ij:'ii> — >■ (^ip^ ijj^jj'^ = 

(^l^l^ + 2n^ $-|-?7i$* and result ((25)) for A which we prove 

shortly, we then obtain the orthogonal HFB equations 
consisting of the modified generalized Gross-Pitaevskii 
equation (GGPE) 



^(r) -^-|-.g(|$|V2n) ) $(r,0 



+gm^*{r,t) - J dr' P{r' ,r,t)<j){r' ,t) 

(15) 

and the orthogonal Bogoliubov-de Gennes equations 
(BdGEs) 



d 

ifi-Q^'^q{r,t) = I dr'L{r,r',t)wq{r',t) 



where we have defined the matrix operator 
L(r,r',t) = 



/:(r,r',t) Miv,v',t) 
M*{r,r',t) -£*{r,r',t) 



and the vector 



yvq{r,t) 



Uq{v,t) 
Vq{r,t) 



(16) 



(17) 



(18) 



for the quasi-particle amplitudes. The operators 
P(r',r,t), £(r, r',t) and A^(r,r',i) are defined by 



P(r',r,i) 
and 



(n(r', r, t)£(r', t) + m(r', r, t)M*{Y', t) 



M{r, r',t) 



Q(r,r',t)£(r',t), 
Q(r, r',t)M{r',t) 



where 

C{r,t) = /io(r)-Ai + 2g(|$(r,0|'+n(r,i) 
M{r,t) EE g{<i>^{r,t) +m{r,t)) , 

with the definitions 

n(r',r,i) = {r^^r' ,t)f^{r,t)) 

= Eg [u*q{r',t)uq{r,t)NBE{eq) 

+vqir',t)v*qir,t) iNBEieq) + l)] 



(19) 
(20) 

(21) 



(22) 
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and 

m(r',r,0 = {f,{r\t)f,{r,t)) 

= Eg [v;ir',t)ug{r,t)NBEieq) (23) 
+ugir',t)v;{r,t) {NBEieg) + l)] 

for the normal and anomalous correlation functions 
n{r',r,t) and m(r',r,i), where 

1 



eigenvalue modes for the orthogonal BdGEs ([27)1 spanned 
by the mode (0, —</>*), let us write wo(r) ~ a(j){r), t'o(r) = 
— /3(/)*(r) for some G C, where <p{r) — ^{r)/y/Nc , 
and 4>{r) is normalized to unity, i.e., 



J dr|(/.(r)|' = l. 



exp(/3e. 



(24) 



is the usual Bose distribution, with /? = l/ksT the tem- 
perature parameter, and ks is the Boltzmann constant. 
Clearly h{r,t) = n(r, r, t) and m{r,t) = m(r, r, t) repre- 
sent the usual thermal and anomalous densities. To prove 
Eq. (fT5|) . we note using results from vector calculus, that 
we can write UM 



|(/dr^t(r,t)LV-(r,t) 



i-/dr{|$|'LyT + |$|'L(|$|2 + 2n) 



-HLVt + 2nL (^|$r + 



(^$*^ -t- m* j Lto - ($2 -I- to) Lto* 
+ (-it + ^jHr,t)Jdr'r){r\t)C*ir',t)r{r',t) 
+ ^r)t (r, t) J dr'r)t (r', i)^ (r' , t)r (r', i)\ L$ 
+ ^ (/dr'ryt(r',t)£(r',i)0(r',t)) r7(r,t) 

+ ^ (/ dr'f7(r', t)>l*(r', t)^(r', t)) r)(r, t)) L$* } 
= lk/rfr(|$|Vn) LVr. 

provided we choose 

/ dr' (f,Hr', t)C{r',t) + f,{r' ,t)M* {r' , t)) (/.(r', t)f,{r, t) 



Then, substituting uq, vq into the first of the modified 
time-independent BdGEs (P7)) . we find that 

eo0 = /dr'Q(r,r') [a£(r',i)</'(r')-/3^(r',t)r(r') 
= a£(r',t)0(r') -/37W(r',t)r(r') 

-0(r) / dr'0*(r) [a£(r', i)0(r') - /37W(r', t)0* (r') 

Multiplying both sides by (p* , and integrating over all 
space, we find [in view of the orthonormality of (jj and of 
the orthogonality condition ^] that 



eo = 0(r)/dr'0*(r) [a£(r',t)(/.(r')-/3>l(r',t)r(r') 
-0(r)/drV*(r) \aC{r' ,t)<j>{r') - l3M{r' ,t)(t)* (r') 
= 0. 



A 



(25) 

thus ensuring that angular and linear momentum are con- 
served. Since standard HFB is a 0-derivable, and hence a 
conserving, theory, particle number and both angular and 
linear momentum conservation hold for standard HFB. 
It can be shown that particle number, and both angular 
and linear momentum conservation also hold for this for- 
malism, the latter being ensured by the inclusion of the 
operator P in (jlSp . The corresponding time-independent 
equations are given by the modified GGPE 



(26) 



We have thus shown in the time-independent case for 
this formalism that there exists a zero-energy eigenvalue. 
However, the existence of a zero-energy eigenvalue does 
not imply that the Pine-Hugenholtz theorem [13|] is sat- 
isfied or indeed that this corresponds to the Goldstone 
mode, nor is the remainder of the energy spectrum "cor- 
rected," and is quite similar to the standard HFB spec- 
trum in spite of the existence of a zero-energy eigenvalue. 

We now obtain from the modified GGPE , the con- 
tinuity equation for the condensate density 

^ mr,t)\^ + V.j(r,t) = n.{rx V) mr,t)\' - ^C{r,t) 
' ot n 

(28) 

where 



i[v,t) = — (<i>(r,t)V$*(r,t) - <i>*(r,t)V$(r,t)) (29) 

is the current density, and where we have defined the 
quantity 



C(r,t) ^ g{rh{v,t)<^*\v,t)-m*{v,t)^^{r,t)^ 
-fG*(r,t) -G(r,t) 



(30) 



+gfh^*{v) - J (ir'F(r',r)0(r') 
and the orthogonal BdGEs 

eqWg{r) = / (ir'L(r,r',)v^fg(r') 



G{r,t) EE ^*{r,t) / dr'P{r',r,mr',t), (31) 



(27) 



where the operators L and P are now time independent. 
To show the existence of a null subspace of zero-energy 



with 



and use this to find an expression for the precessional fre- 
quencies of off-axis vortices and vortex arrays in quasi- 
two-dimensional BECs. To do this we consider a station- 
ary vortex system in the frame rotating at the preces- 
sional frequency of the vortex or vortices. 
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III. CALCULATION OF PRECESSIONAL 
FREQUENCIES OF VORTICES IN 
QUASI-TWO-DIMENSIONAL BEGS 

In the quasi-2D regime, the time-indcpcndcnt orthog- 
onal HFB equations in a frame rotating with angular 
frequency O are given in polar coordinates by the time- 
independent GGPE 



2D 



2n $ + C2Dm^ 



j^r'dr'de'P{r',e',r,e)^{r',e') 



(32) 



and the 2D time-independent Bogoliubov de Gennes 
Equations (BdGEs) 



where 
and 

t{r,ey,6') 



lo^ /o" r'dr'de't{r, 9, r' , 0')wg(r', 0') 



Uq{r, 9) 
vg{r,9) 



(34) 



C{r,9,r',9') M{r,9,r',9') 
-M*{r,9,r',9') -t {r,9,r' ,9') 



Here we have defined the operators 

C{r,e,r',e') = Q{r,e,r',e')t{r',e') 
M{r,9,r',e') = Q{r,e,r' ,e')M{r' ,9') 

and 



(35) 
(36) 



P{r',9',r,9) = {h{r' ,9' ,r,9)C{r' ,9') 

+m{r',9',r,9)M*{r',9')) /^/NJ{) 

(37) 

with 

C{r,9) = hn{r,9)- fi + 2C2D[mr,9)f +n{r,9)^ 
Mir,9) = C2D{^^{r,9)+m{r,9)), 

(38) 



Q{r, 9, r', 9') = -S{r - r')S{9 - 9') - 4>{r, 9)^*{r', 9'), 

(39) 

and 

f^nir,9) = -(^ + l§-^ + ^^)+ial^W. (40) 
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We note that in the quasi-2D regime that SI reduces to 
the component = fl^ along the 2;-axis. We use the 



dimensionless units ''o = y for the length scale, 

and to = for the time scale, and express the con- 
densate w ave fu nction and quasiparticle amplitudes in 
units of -y/ N/ rg . Then the energies arc given in har- 
monic oscillator units hoJr/2. The nonlinearity constant 



C2D is found by expressing the interaction term g in 
dimensionless units, and integrating over all z yielding 

1 /2 

C2D = Stt (^) where A'' is the number of atoms, 

Us is the scattering length, and A = uJz/oJr is the trap 
aspect ratio. 

Let R{r) = Rc($(r)) and /(r) = Im($(r)). Then in 
the two-dimensional system, the continuity equation for 
the condensate density is given by 



nr^A{r,9) = B{r,9) 
where we have defined 

Air,9)^[R- + I- 

and 



(41) 



(42) 



B{r,9) = 



a' I 



+ [R-902 - I 



09^ 



)]+i^C(r,^). 



with 

C{r,9) 



(43) 



C2D (rhir,9)<i>*\r,9) -7fi*{r,9)^^ir,9) 
+G*{r,9)-G{r,9) 



and 

G{r,9) = $*(r,( 



(44) 



27r /"Cso 



r'dr'd9'P{r',9',r,9)(j>{r',9'). 

(45) 

Multiplying both sides by A{r,9), and integrating over 
all space gives us the expression 



JT ^(^' ^)-^(^' 0ydrd9 
Cj,^r^A{r,9frdrd9 



(46) 



which we solve self-consistently with the time- 
independent orthogonal HFB equations in the frame ro- 
tating at angular frequency Q. The converged solu- 
tion then represents a stationary solution in the rotating 
frame, and £7 corresponds to the precessional frequency 
of the off-axis vortex or vortex array. 

A vortex at position (ri, ^1) may be created by expand- 
ing the condensate wave frmction in terms of modified 

Laguerrc basis functions |x;^^(f', ^)| as follows: 
Define 



5W = {Z,n| (I,n)e5-{(1,0)}} 
where S = {l,n \ n = 0, . . . ,1 = 0, ±1, . . .}, then 
^r,9)= J2 ^inx\l\r,9) 



(47) 



(48) 



provided (ri, ^1) is not a root of ^inir, 9), where we define 

c;w{ri,tii) 



5 



with 



0.7 



iW 



2nl 



2^ V(^+I^l)! 



1/2 



-^'/'r^LWir^ (50) 



and in (a;) a modified Laguerre polynomial of order n. 
The superscript in the representation of the modified ba- 
sis function, indicates a single vortex. This is motivated 
simply by expanding the condensate wave function in 
terms of the complete Laguerre basis {^in{r, 0) \ l,n G S} 
and considering that $(^,6*1) = Y^inain^i„{ri,0i) = 0. 
Imposing this condition implies dependency of one of the 
coefficients ain on the others, reducing the dimension of 
the vector space by one. Here we choose to represent aio 
in terms of the remainder {a;„ \ l,n€ S^^^}. We find 
that aio = -J2i,nesw {^in{ri,di)/(,io{ri,9i)) ain, hence 
we can write $(r, 9) in the form (|48|) . 

This procedure may be extended to vortices 
{{ri,9i) , . . . , (''Ar„, ^Af„)} by an iterative process, writing 



$(r, ( 



(51) 



where we have defined 



x\:!^\r,9) . xif-^^(.,^)- 4;:!;;^"-'"-; x!vto^(^,^). 

Xn^.o (rN^,9Nj 

and where we have generalized the definition (|47p in the 
single- vortex case to 

={l,n\ {I, n)eS- {(1, 0) , . . . , (iV„, 0)}} (53) 

in the multi-vortex case. 

We illustrate the method by considering a dilute Bose 
gas consisting of 2000 *^Rb atoms in an axially symmet- 
ric harmonic trap, with radial and axial harmonic trap- 
ping frequencies of fir = 2ttx lOHz and ft^ = 2ttx 400Hz, 
respectively. Thus the axial confinement is sufficiently 
strong that all excited axial states may be neglected, 
hence the axial (z) component can then be integrated 
out, and the BEC may effectively be treated as two- 
dimensional. In order for this quasi-two-dimensional 
approximation to hold, we require that the axial har- 
monic transition energy huJz be much greater than the 
temperature T and the mean-field energy of the BEC 
Here the axial harmonic energy transition corre- 
sponds to a temperature of Ta^ — hojz/kB — 19.2nK, 
and the mean-field energy to an effective temperature of 

Tmf = CaD^'^Vfcs = •84nK, where ks is Boltzmann's 
constant, and fc^ = 2kB /fujJr- In all these simulations the 
BEC temperature T < lOnK, so this approximation is 
justified. 

Stationary solutions are found in the case of a sin- 
gle off-axis vortex, and in the cases of triangular and 
hexagonal vortex arrays, provided the interactions be- 
tween the votices are not too strong. In cases where 
there are strong interactions between vortices, approx- 
imately stationary solutions can be found, establishing 



(a) 



(b) 



♦ ♦ 




(d) 



3.5 



T(nK) 



FIG. 1: Triangular vortex array, (a) Precessional frequency 
Q versus lattice parameter a at temperatures 1, 2, 5 and 
7.5 nK (triangles, diamonds, circles, squares); (b) preces- 
sional frequency fl versus temperature T for lattice parameter 
1.5\/3, 1.65^, 1.8^, and 2^/3 (t riangles, diamonds, circles, 
squares). Hexagonal vortex array, (c) Precessional frequency 
Q (in units of radial trapping frequency oj, ) versus lattice 
parameter a at temperatures 1, 2, 5 and 7.5 nK (triangles, 
diamonds, circles, squares); (d) precessional frequency fl ver- 
sus temperature T for lattice parameter 2.75, 2.8, 2.9 and 3 
(triangles, diamonds, circles, squares). All frequencies are in 
units of the radial trapping frequency ujr. All positions are in 
units of the harmonic oscillator length ro. 



lower and upper bounds for the precessional frequency, 
and the precessional frequency determined more accu- 
rately in time-dependent simulations. 

Precessional frequencies as a function of vortex posi- 
tion a and of temperature T for the single off-axis case 
are presented in Figs. 1(c) and (d) of a previous publi- 
cation [1^, and will not be repeated here. Instead, we 
concentrate on the cases of triangular and hexagonal vor- 
tex arrays. The precessional frequency i7 as a function 
of temperature T for various values of lattice parameter 
a and of lattice parameter a at various temperatures are 
shown in Figs. 1(a) and 1(b) for the triangular vortex 
array and in Figs. 1(c) and 1(d) for the hexagonal vor- 
tex array. We see that decreases with increasing lattice 
parameter a over the range shown (within the Thomas- 
Fermi radius). This is in accord with the areal density 
law given by the approximation Q 



m 



1 



r?2 



2vr [Rl^ - r^Y 



In 



(54) 



where is the areal density of vortices approximated 
by n„ — Ny/{TrR'^p), is the number of vortices in 
the vortex array, and Rtf is the Thomas-Fermi radius 
and may be understood by considering that the vortex 
precession velocity is determined primarily by the back- 



ground velocity field around the vortex core, and the core 
shape [H, [13] • The effect of their image charges in the 
condensate boundary is less important, in contrast with 
the single vortex case [3) [3 ■ In the single off-axis case 
il increases with increasing a p^ . This is in agreement 
with GPE results [3, [l^ and may be explained in terms 
of the effect of the vortex's image charge in the boundary 
of the condensate 18, 19]. 

Figures 2(a) and 2(b) show, respectively, simulated ab- 
sorption ima ges for the condensate density and the ther- 
mal density [20|, and Fig. 2(c) the condensate phase 
in the xy-plane for a single off-axis vortex situated at 
a = 0.5. We note the phase singularity in Fig. 2(c) in- 
dicating a singly charged vortex. Similarly Figs. 2(d) 
and 2(e) show, respectively, simulated absorption images 
for the condensate density and the thermal density, and 
Fig. 2(f) shows the condensate phase in the xy-plane 
for a triangular vortex latice having lattice parameter 
a = 1.65a/3. Figs. 2(g)-2(i) show the respective plots 
for a hexagonal vortex latice having lattice parameter 
a = 2.9. The phase singularities in Figs. 2(f) and 2(i) 
indicate respectively, three and seven singly charged vor- 
tices. We note from the simulated absorption images for 
the condensate in Figs. 2(a), 2(d) and 2(g) that the con- 
densate extends further outward for the triangular and 
hexagonal vortex arrays (i.e., larger Thomas-Fermi ra- 
dius) than for the single vortex case. This is consistent 
with our findings concerning the dependency of on a, 
where we find that decreases with increasing a in the 
vortex array case, in contrast to the single off-axis vortex 
case. Unfortunately, it is not possible to find stationary 
solutions in the limit of sufficiently large numbers of vor- 
tices to validate against the areal density approximation 

0]; nevertheless we perform this calculation for the 
case of seven vortices. 

Figures 3(a) and 3(b) represent plots of processional 
frequency versus lattice parameter a at temperatures 
T = 2 nK and T = 5 nK, respectively, for vortex arrays 
composed of two, three, and seven vortices. We see that 
fl increases with an increasing number of vortices. The 
inserts of Figs. 3(a) and 3(b) seem to indicate that 51 
tends asymptotically to an upper limit as the number of 
vortices increases, presumably the limit predicted by the 
areal density approximation (1541) . In the dimensionless 
units used here, this may be written 



1 



R2 

IXrp p 



1 



■In 



(55) 



Therefore, for the hexagonal array where Ny = 7, and 
taking Rtf ~ 3.5ro, we predict to first approximation a 
precessional frequency of £7 ~ .57wr- Taking into ac- 
count the second term in (|55p with r ~ 0.5Rtfi we 
predict 51 ^ .550;^ (so the second term is unimportant 
in the regime considered here). This is of the order of 
20% lower than the value predicted using the continuity 
equation calculation. We note that the areal density ap- 
proximation is valid only for a large number of vortices. 
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FIG. 2: Simulated absorption images ^0] for ofT-axis vortex, 
triangular vortex array and hexagonal vortex array showing 
(a) condensate density, (b) thermal density, (c) condensate 
phase for a single vortex at a = 0.5, and (d) condensate den- 
sity, (e) thermal density, (f) condensate phase for a triangular 
vortex array with lattice parameter a = 1.65\/3, and (g) con- 
densate density, (h) thermal density, (i) condensate phase for 
a hexagonal vortex array with lattice parameter a —2.9. All 
positions are in units of the harmonic oscillator length ro. 
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FIG. 3: Precessional frequency Q versus lattice parameter a 
for vortex arrays composed of two (circles), three (triangles), 
and seven vortices (squares) with cubic spline interpolation 
(dashed line) (a) at T = 2 nK, and (b) at T = 5 nK. The insert 
in each respective figure shows the precessional frequency 
versus number of vortices for vortex arrays composed of two, 
three, and seven vortices having lattice parameter a = 3.0 
with cubic spline interpolation (dashed line). All frequencies 
are in units of the radial trapping frequency tOr. All positions 
are in units of the harmonic oscillator length ro. 



Furthermore, referring to the inserts of Figs. 3(a) and 
(b) reveals that the precessional frequency approaches 
an asymptotic value of ^ .72. We conclude, however, 
that these results are consistent in the asymptotic limit. 
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IV. EVOLUTION OF TIME-INDEPENDENT 
SOLUTIONS 

Time-dependent simulations using the time- 
independent data initial state are in good agreement 
with the above calculations. For a single vortex at 
a = 0.5 radial harmonic oscillator units from the axis, 
the processional frequency was estimated using a least- 
squares fit of a sinusoid to the vortex x-displacement 
versus time and was found to he Qls ~ 0.3794a;r, in 
good agreement with the value of 17 = 0. 3727 uJr, as pre- 
dicted in the time-independent calculations. The error 
in the prediction scales as the number of computational 
basis functions which, for practical reasons, is only 
209 in these calculations. For 839 computational basis 
functions, the time- independent calculations predict a 
value of = 0.3761a;r- We see no evidence of dissipation 
during the time of the simulation of five trap cycles, 
and the results reveal that the trajectory of the vortex 
constitutes a circle. 

Time-dependent simulations were also performed at 
T = 5 nK for a symmetrical triangular vortex array with 
lattice parameter a = 1.65\/3 (i.e., three vortices sym- 
metrically positioned at 1.65 radial harmonic oscillator 
units from the axis), and for a symmetrical hexagonal 
triangular vortex centered on the axis having lattice pa- 
rameter a = 2.85 (i.e., seven vortices), again using time- 
independent calculations for the initial state. In both 
simulations, the off-axis vortices precess in a circle, and 
again we see no evidence of dissipation during the time 
of the simulation of five trap cycles. The precessional fre- 
quency fl^s = 0.5856a;r for the triangular vortex array 
was estimated using a least squares fit as before, in good 
agreement with the value of = 0.59380;^, as predicted 
in the time-independent calculations for 209 computa- 
tional basis functions. 



V. STIRRING OF THE CONDENSATE 

In these, and in other time-dependent simulations, 
the temperature T remains fixed. It is not necessary 
to vary ^ to maintain particle conservation since the 
rate of change in the condensate population dNc/dt and 
in the non-condensate population dN/dt are equal and 
opposite, depending only on m and $, i.e., dNc/dt — 

—dN/dt — i| / dr — m*<I>^^ , therefore particle 

conservation always holds. In these simulations we use 
one or more Gaussian rotating potentials in order to stir 
the condensate. The stirring imparts angular momen- 
tum to the condensate, and therefore one would expect 
the production of vortices in cases where the angular fre- 
quency of the stirrer(s) exceeds a certain critical value. 
This should correspond (at least approximately) to the 
velocity of the stirrer with respect to the fluid exceeding 
the local Landau critical velocity (i.e. the local speed of 
sound). In what follows we shall establish a critical stir- 



ring frequency, below which no vortices are produced in 
regions of appreciable density (i.e., within the Thomas- 
Fermi radius). In these simulations we find that this 
critical frequency corresponds to a velocity of the stirrer 
through the BEC that is slightly in excess of the local 
Landau critical velocity (see Table 1). In the simulations 
where the critical stirring frequency is exceeded, we are 
able to determine by the method of least squares, as al- 
luded to above, the precessional frequencies and compare 
these with the values predicted in the time-independent 
calculations using the continuity equation for the con- 
densate density. We shall see presently that the agree- 
ment between the time-dependent values and the pre- 
dicted values using the time-independent calculations is 
good, as one would expect in the adiabatic case since 
then the term ih-^ \^\^ in the time-dependent continuity 
equation for the condensate density (p8)) would be small, 
i.e., of order of the fluctuations in the chemical poten- 
tial. We stir the BEC using a single Gaussian stirrer 
at a— 1.5ro off axis, of amplitude 10 hur/'i, full width 
at half maximum (FWHM) 0.82 ro- Two cases are con- 
sidered. In the flrst we stir the BEC at a sub-critical 
stirring frequency 0.38 tUr, with the stirrer on for 10 trap 
cycles. In the second case we stir the BEC with stirring 
frequency Q — 0.5^^ , with the stirrer on for 4.5 trap 
cycles. In both cases the stirring frequency is ramped 
up from zero to the specifled value over 0.2 trap cycles. 
In the second case two vortices are produced, and the 
simulated absorption images are shown in Fig. 4, and 
the trajectories of these vortices in Fig. 5 |21|]. Figure 
6 1(a) shows the condensate, thermal, and total angular 
momentum in the first case where the stirring frequency 
(ri = 0.38w,.) is sub-critical, and Fig. 611(a) in the sec- 
ond case where = 0.5LOr which is above the critical 
stirring frequency, and where two vortices are produced 
as shown in Figs. 4 and 5. In the case of the first sim- 
ulation, the rotational frequency is sub-critical, so no 
vortices are produced within the Thomas Fermi radius, 
therefore negligible angular momentum is transferred to 
the condensate. This angular momentum is due to cir- 
culation corresponding to vortices in regions of negiible 
density. We see in both cases that a small amount of an- 
gular momentum is also transferred to the thermal cloud. 
We note that remains fixed once the stirring ceases 
after 10.2 trap cycles in the first case, and 4.7 trap cy- 
cles in the second in accordance with the conservation of 
angular momentum. We note in the first case where the 
stirring frequency is subcritical that the angular momen- 
tum for the thermal population increases for the first few 
trap cycles, reaching a saturation value after ~ 4 — > 5 
trap cycles. Simulations over a much longer time scale 
(not presented here) reveal a slow cyclical change in the 
thermal population angular momentum of a few percent- 
age points. Once the stirrer is switched off, we see a 
steady increase in the thermal population angular mo- 
mentum with a resultant loss in the angular momentum 
of the condensate (since the overall angular momentum 
is conserved). During the time of stirring, the condensate 
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FIG. 4: Simulated absorption images ^20|1; stirring of con- 
densate in an counterclockwise direction using one Gaussian 
stirrer at a= 1.5 ro off axis, of amplitude 10 hjjJr/2, FWHM 
0.82 ro, switched on adiabatically over 0.2 trap cycles, stirring 
frequency 0.5 Ur, stirrer on for 4.5 trap cycles after (a) 0.5, 
(b) 1, (c) 1.5, (d) 2, (e) 2.5, (f) 3, (g) 3.5, (h) 4, (i) 4.5, (j) 5, 
(k) 5.5, (1) 6, (m) 6.5, (n) 7, (o) 7.5, (p) 8, (q) 8.5, (r) 9, (s) 
9.5, and (t) 10 trap cycles. All positions are in units of the 
harmonic oscillator length ro. 



angular momentum varies cyclically corresponding to the 
variation of the distance of vortices from the axis (see also 
Ref. t22j). In the case where the stirring frequency is 
subcritical, these vortices would lie outside the Thomas- 
Fermi radius, and would therefore be subthreshold. Figs. 
61(b) and 611(b) show, respectively, the changes in con- 
densate and thermal populations, where we see rigorous 
observance of particle number conservation. In order 
to understand why there exists a critical frequency for 
the production of vortices within the Thomas-Fermi ra- 
dius, we calculate the Landau critical velocity v^c given 
by the local speed of sound at the point of the stirrer 
s = -y/ ng/m, where n is the density of the condensate 
and m is the mass of a particle. We use the local density 
approximation and calculate the speed of sound s at the 
stirrer using the local density n of the superfluid at the 
point of the stirrer. In these simulations, we measure the 
two-dimensional density n*^^^^ where we have integrated 
out the axial component, assuming the axial confinement 
is sufficiently tight that all axial modes except the low- 
est can be ignored, and we can write the local Landau 
critical frequency in terms of the local two-dimensional 
density n^^^' as Vlc = VsTScW*^- The speed of 
the stirrer is given by Ustir = ^ttTs^s for a stirrer at ra- 
dius rotating at angular frequency fis. We can rewrite 
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FIG. 5: Stirring of condensate using one Gaussian stirrer 
at a= 1.5ro off axis, of amplitude 10fej,./2, FWHM 0.82ro, 
switched on adiabatically over 0.2 trap cycles, stirring fre- 
quency Q.^u)r, stirrer on for 4.5 trap cycles (a) x-displacement 
of vortex 1, (b) y-displacement of vortex 1, and (c) trajectory 
of vortex 1 over a period of 10 trap cycles (d) i- displacement 
of vortex 2, (e) ^-displacement of vortex 2, and (f) trajectory 
of vortex 2 over a period of 10 trap cycles. All positions are 
in units of the harmonic oscillator length ro. 





Vstiv/ Vlc 


Lz(iVfi) 


0.38 


1.12 


0.035 


0.39 


1.15 


0.045 


0.4 


1.18 


0.5 (varies between ~ 0.25 and ~ 0.7) 


0.45 


1.33 


0.55 (varies between ~ 0.3 and ~ 0.8) 


0.5 


1.48 


0.6 for one vortex, ~ 0.8 for two vortices 
and ~ 2.2 for three vortices 



TABLE L Relationship of the local stirrer velocity with the 
average z-component of the angular momentum Lz for time- 
dependent HFB simulations. The Landau critical velocity is 
that corresponding to the density at the stirrer location, and 



,(2D) 



0.02 r 







this in terms of the Landau critical velocity as 

-.VLC- (56) 



Table 1 shows numerically calculated values of the time- 
averaged 0-component quantum expectation value of the 
angular momentum Lz and how this is related to the 
local stirring velocity in units of the Landau critical ve- 
locity for the stirring angular frequencies 0.39, 0.4, 0.45, 
and 0.5 Ur- We see evidence of a critical stirring angular 
velocity between 0.39 and 0.40;^ corresponding to a lo- 
cal stirrer velocity just slightly above the Landau critical 
velocity. These results are consistent with simple GPE 
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FIG. 6: (Color online) Stirring of condensate (I) using one 
Gaussian stirrer at a= 1.5ro off axis, of amplitude 10 hu}r/2, 
FWflM 0.82 ro, switched on adiabatically over 0.2 trap cy- 
cles, stirring frequency 0.38 LUr, stirrer on for 10 trap cycles, 
and (II) using one Gaussian stirrer at a= 1.5ro off axis, of 
amplitude 10 hLJr/2, FWHM 0.82 ro, switched on adiabati- 
cally over 0.2 trap cycles, stirring frequency 0.5 ujr, stirrer on 
for 4.5 trap cycles showing (a) z component of total angu- 
lar momentum (solid line), condensate angular momentum 
(dashed line), and thermal population angular momentum 
(dash-dotted line) versus time, (b) change in thermal pop- 
ulation and condensate population (solid line) versus time. 
The change in the populations AA'^ is given in units of the 
number of atoms A'^. All angular momenta are in units of Nh. 



simulations (results not presented here). 

No experiments have been performed on stirring a BEC 
using a narrow rotating probe as considered here but ex- 
periments [23] have been performed that show evidence 
of a critical stirring angular velocity using elliptical de- 
formation of the trap. Experiments to test the superfluid 
critical velocity using the linear motion of a narrow probe 
laser beam have also been performed (23 - [26| , finding crit- 
ical velocities for the onset of dissipation ~ 5 — 10 times 
smaller than the critical velocity inferred from the bulk 
Bogoliubov speed of sound. Time-dependent GPE sim- 
ulations ^] also indicate that dissipation first occurs at 
speeds significantly less than the speed of sound. How- 
ever, the linear motion of the probe in these works excite 
vortex dipole or solitons, and are not therefore directly 
applicable to the present work. 

In a previous study p^ . similar stirring simulations 
were carried out using time-dependent GPE simulations. 
However, the critical stirring frequency referred to in 
these results are the thermodynamic critical frequency 



(at which the vortex-free and central vortex states have 
equal energies), whereas the results presented here re- 
fer to the critical frequency at which vortices are cre- 
ated within the Thomas-Fermi radius. Nevertheless, the 
results of Rcf. [22] are consistent with the simulations 
performed in this work. 

The precessional frequencies for single off-axis vortices 
were obtained for simulations deploying the single stirrer 
atn = 0.5LOr for 2.5 trap cycles. Least squares analysis of 
the vortex trajectories as described above reveals that the 
vortex in the second simulation developed at a distance 
of ^ 1.7 harmonic oscillator units from the axis, with a 
precessional frequency of ^ 0A2u!r, in very good agree- 
ment with the off-axis time-independent calculations (see 
Figs. 1(c) and 1(d) in Ref. [11]). 

In order to test the predicted precessional frequencies 
for triangular vortex arrays, three equi-spaced Gaussian 
stirrers at a= 1.5 rg off axis, of amplitude 10 hujr/2 
(i.e. of the same order as the chemical potential, /i > 
lOhojr/^), FWHM 0.82 rg, stirring frequency ramped up 
from zero to 0.50;^. over 0.2 trap cycles, are used to stir the 
BEC for 2.5 trap cycles. Thus a symmetrical triangular 
vortex array was produced. The least-squares estimates 
of the precessional frequency of f^LS = O.TuJr, are in very 
good agreement with the interpolated values from Figs. 
1(a) and 1(b) of r2 - 0.7 ujr. 

The above results for the stirring simulations indicate 
very good agreement with the values predicted using the 
continuity equation for the condensate density in the 
time-independent calculations. These results were at- 
tained in the case where the stirrers were introduced adi- 
abatically. It should be noted that the reason for using 
three Gaussian stirrers, as opposed to just one, is not 
that we need three stirrers to create three vortices (one 
stirrer is sufficient for this purpose provided the BEC is 
stirred long enough at a sufficiently high frequency) , but 
rather in order to create a symmetrical triangular vortex 
array so the time-independent predictions can be tested. 



VI. BREAKING OF THE AXIAL SYMMETRY 
OF THE TRAPPING POTENTIAL 



Noether's theorem |28| states that in a conservative 
physical system, a differentiable symmetry of the La- 
grangian (and, it can be shown, of the Hamiltonian) has 
a conserved quantity associated with the system (i.e a 
corresponding conservation law). In this case the axial 
symmetry of the trapping potential implies the conserva- 
tion of the axial component (i.e., the z-component) of the 
angular momentum. The quantum expectation value of 

the angular momentum is given by (L) = J dr ^i/j^Lt/;^, 
where ijj — ^ + fj as above. Thus using the expression 
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time (trap cycles) 

FIG. 7: (Color online) Breaking of axial symmetry; (I) for 
single off-axis vortex for turn-on time of 20 trap cycles for ec- 
centricity parameter e= 0.25, and (II) for single off-axis vortex 
for turn-on time of 20 trap cycles for eccentricity parameter 
e= 0.1 showing (a) z-component of total angular momentum 
and condensate angular momentum (lower line) versus time 
and (b) change in thermal population and condensate popu- 
lation (lower line) versus time. The change in the populations 
AA^ is given in units of the number of atoms A'^. All angular 
momenta are in units of Nti. 



(|22p for the non-condensate density matrix, we find that 

(L) = -ih Jdr $* (r X V) $ + J2q ("g", (r x V) 
+ (71, + 1)1,, irxV)v;)]. 

(57) 

It can be shown from the time-dependent HFB equations 
TB and SB that 



dt 



= / dr 1$ 



LT/t. 



(58) 



This determines the rate of change in angular momen- 
tum in the BEC system, and we see that -^L^ — for an 
axially symmetric trapping potential in accordance with 
Noether's theorem. The implication of this is that when 
this symmetry is broken, Lz is no longer conserved, and 
since there is no external source of angular momentum 
here, we conclude that angular momentum will be lost 
since jf.Lz will now be non-zero by equation This 
leads to vortex decay, and the vortex spirals outward as 
angular momentum is lost. This is in accordance with 
work done by Zhuravlev et. al. [1^ on the arrest of 
vortex arrays by a trap anisotropy, although this work 
deals with large vortex arrays, whereas the the calcula- 
tions performed here concern single off-axis vortices and 
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FIG. 8: Simulated absorption images teu] for breaking of axial 
symmetry with turn-on time of 20 cycles for single off-axis 
vortex for eccentricity parameter e= 0.25 after (a) 8, (b) 16, 
(c) 24, (d) 32, (e) 40, (f) 48, (g) 56, (h) 64, (i) 72, (j) 80, (k) 
88, (1) 96, (m) 104, (n) 112, (o) 120, (p) 128, (q) 136, (r) 144, 
(s) 152, and (t) 160 trap cycles. All positions are in units of 
the harmonic oscillator length ro. 



small vortex arrays. In order to test this, we introduce an 
eccentricity into the trapping potential, thereby breaking 
the axial symmetry. Initially the trapping potential is ax- 
ially symmetric and the eccentricity introduced by means 
of the perturbation potential 



— cos' 



e) 



(59) 



which is introduced adiabatically after two trap cycles 
over a period of 20 trap cycles. In these simulations, we 
introduce eccentricities of e = 0.25 and e — 0.1 in both 
the cases of a single off-axis vortex and of a triangular 
vortex array. In the single vortex case, the initial state 
consists of a BEC with a single precessing vortex situ- 
ated at position a — 1.1 trap units from the axis. The 
z-component of the angular momentum versus time 
in trap cycles is plotted in Figs. 71(a) and 711(a) for the 
cases of e = 0.25 and e = 0.1, respectively, where we see 
a decrease in Lz with time once the symmetry of the 
trapping potential has been broken. Figures 71(b) and 
711(b) show the respective changes in the thermal and 
condensate populations, and reveal an increase in the 
thermal population as the vortex spirals outwards and 
finally decays. This is consistent with the notion that 
the vortex decays into excitations when it reaches the 
condensate boundary [s^l . Simulated absorption images 
for e = 0.25 with turn-on time of 20 trap cycles are shown 
in Figs. 8(a)-8(t). In the simulations with e = 0.25 the 
lifetime of the vortex is approximately 145 to 150 trap 
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FIG. 9: (Color online) Breaking of axial symmetry, (I) for 
triangular vortex array for turn-on time of 20 trap cycles for 
eccentricity parameter e= 0.25, and (II) for triangular vor- 
tex array for turn-on time of 20 trap cycles for eccentricity 
parameter e= 0.1 showing (a) z-component of total angular 
momentum and condensate angular momentum (lower line) 
versus time and (b) change in thermal population and change 
in condensate population (lower line) versus time. The change 
in the populations AA^ is given in units of the number of 
atoms A''. All angular momenta are in units of Nh. 



cycles, and for e = .1 substantially larger, with an ex- 
trapolated value of the order of 700 trap cycles. The vor- 
tex precessional frequencies are in good agreement with 
the time-independent calculations, and we note that in 
all cases from the vortex trajectories that the vortices 
precess faster as they move outward [31], which, again, 
is in agreement with the time-independent results (see 
Fig. 1(c) of Ref. [HI), before finally decaying. This is 
also consistent with the reduction in angular momentum 
[see Figs. 71(a) and 711(a) for the cases of e = 0.25 and 
e = 0.1 respectively]. The vortex moves outwards and 
precesses faster as it does so in agreement with Fig. 1(c) 
of Ref. [l^ , accompanied by a loss in angular momentum 
[see Figs. 71(a) and 711(a)]. As discussed above, this loss 
of angular momentum is due to the symmetry breaking 
of the trap and the loss of the angular momentum is as 
given by (I55)) . The loss of the vortex after ~ 145 — >■ 150 
trap cycles (loss of tracking after ~ 146 trap cycles) for 
the case e = 0.25 is consistent with the decline in Lz seen 
in Fig. 7 1(a) and the sudden increase in the thermal pop- 
ulation in Fig. 7 1(b) at ~ 150 trap cycles. These results 
are also in qualitative agreement with recent simulations 
using a classical field treatment comprising the projected 
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FIG. 10: Simulated absorption images [20|] for breaking of ax- 
ial symmetry for triangular vortex array over 20 trap periods 
with lattice parameter a = 1.65\/3 for eccentricity parameter 
e= 0.25 after (a) 2.5, (b) 5, (c) 7.5, (d) 10, (e) 12.5, (f) 15, 
(g) 17.5, (h) 20, (i) 22.5, (J) 25, (k) 27.5, (1) 30, (m) 32.5, (n) 
35, (o) 37.5, (p) 40, (q) 42.5, (r) 45, (s) 47.5, and (t) 50 trap 
cycles. All positions are in units of the harmonic oscillator 
length To. 
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FIG. 11: Breaking of axial symmetry for triangular vortex 
array for turn-on time of 20 trap cycles for eccentricity pa- 
rameter e = 0.25 showing (a) x displacement of vortex 1, (b) 
y displacement of vortex 1, (c) trajectory of vortex 1, (d) x 
displacement of vortex 2, (e) y displacement of vortex 2, (f) 
trajectory of vortex 2, and (g) x displacement of vortex 3, (h) 
y displacement of vortex 3, and (i) trajectory of vortex 3. All 
positions are in units of the harmonic oscillator length tq. 
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Gross-Pitaevskii equation [32] ■ 

Similar simulations were performed for the triangular 
vortex array with lattice parameter a — 1.65-\/3j where 
the trapping potential is initially axially symmetric, and 
the eccentricity is ramped up to e = .25 and e = .1 over 
a period of 20 trap cycles. 

The angular momentum calculated using Eq. (j57p is 
plotted in Fig. 9 1(a) for the case of e = 0.25 and in Fig. 
9 11(a) for the case of e = 0.1. Figures 91(b) , and 911(b) 
show the respective changes in the thermal and conden- 
sate populations and reveal increases in thermal popu- 
lation as each successive vortex spirals outward and fi- 
nally decays, again in accordance with decay of the vortex 
into excitations when it reaches the condensate bound- 
ary [3^. The simulated absorption images for e = .25 
are shown in Figs. 10(a)-10(t), and the vortex trajecto- 
ries in Figs. ll(a)-ll(i). We note from the trajectory of 
vortex 3, shown in Figs. 11(g) and 11(h), the decrease 
in the precessional frequency following the departure of 
a vortex, in agreement with the time- independent pre- 
dictions for vortex arrays shown in Figs. 3(a) and 3(b), 
which show the dependence of the precessional frequen- 
cies of the vortex array on the lattice parameter for the 
cases of three vortices (our initial configuration) , and of 
two vortices, and also with Figs. 1(a) and (b) which give 
the time-independent predicted precessional frequencies 
for the triangular vortex array at various temperatures 
and values of lattice parameter. The prediction here is 
n ~ .6uJr for lattice parameter a ~ 2-\/3 — 3.46 and is in 
good agreement with the measured time-dependent value 
[see Figs. 11(a) and 11(b) for Xpos versus T and ypos ver- 
sus T and the time- independent predictions in Figs. 1(a) 
and (b)]. We also see good agreement with the measured 
value of the precessional frequencies of a single vortex (af- 
ter the other two vortices have left the condensate) indi- 
cated in Figs. 11(g) and 11(h), and the time-independent 
predictions given in Figs. 1(c) and 1(d) in Ref. [l5|. We 
note that Figs. 10 and 11, which show the vortex de- 
cay in the case e = .25, are consistent indicating a decay 
time of ~ 24 trap cycles for the first vortex and a de- 
cay time of ^ 38 trap cycles for the second. This is also 
consistent with Fig. 91 which shows sudden changes in 
the thermal population at these times, and a steady de- 
crease in Lz with slight drops accompanying the decay of 
each vortex. We also note a sudden change in the ther- 
mal population after 30 trap cycles. This accompanies 
the complex precessional motion of the third vortex at 
^ 30 trap cycles as shown in Figs. 11(g) and 11(h) due 
to the interaction between the vortices. The lifetimes 
of the vortices in the case of the triangular vortex for 
e = .25 vary considerably - the first vortex to leave the 
condensate has a lifetime of between ~ 20 and ~ 25 trap 
cycles, while the second has a lifetime varying between 
^ 25 and ^ 38 depending on the scenario. In all cases 



the loss of a vortex is characterized by a slight drop in 
Lz and sudden increases in the thermal population with 
a corresponding drop in the condensate population (see 
Fig. 9). The departure of the vortices is also marked by 
the trajectories of the vortices shown in Figs. ll(a)-ll(f) 
where we see the end of the vortex tracks. The lifetime 
of the remaining vortex can be inferred from Fig. 9 1(a) 
only by extrapolation and is of the order of 300-400 trap 
cycles [since the residual angular momentum (i.e. the 
angular momentum at the time that only a single vortex 
remains) is higher than in the single vortex case] . In the 
case of the triangular vortex array for e = 0.1, the life- 
time of the first vortex is of the same order (^ 25 — 30 
trap cycles), but the lifetimes for the second and third 
vortices are considerably longer, with the second vortex 
having an extrapolated lifetime of ^ 400 trap cycles [see 
Fig. 9 11(a)]. It is impossible to ascertain the lifetime 
of the remaining vortex from the simulations, but it is 
probably of the order of a few thousand trap cycles that 
is impractical to simulate with the facilities available. 



VII. SUMMARY 

Wc have developed a novel method utilizing the 
continuity equation for the condensate density to make 
predictions of the precessional frequency of single 
off-axis vortices, and of vortex arrays in Bose-Einstein 
condensates (BECs) at finite temperature. We also 
presented an orthogonalized Hartree-Fock-Bogoliubov 
(HFB) formalism for which a zero-energy eigenvalue 
exists, in contrast to the standard HFB theory. We 
solved the continuity equation for the condensate density 
and the time-independent orthogonalized HFB equations 
self-consistently in the frame rotating at the preces- 
sional frequency. Thus we were able to find stationary 
solutions for quasi-two-dimensional rotating systems 
in the corotating frame. We compared these results 
with time-dependent predictions where we simulated 
stirring of the condensate, finding good agreement with 
the predicted precessional frequencies. These results 
are completely consistent with GPE simulations, and 
are in qualitative agreement with the projected Gross- 
Pitaevskii simulations and with the areal density law. 
We also verified that angular momentum is conserved 
in the quasi-two-dimensional system for an axially 
symmetric trapping potential and that breaking this 
symmetry leads to the loss of angular momentum and 
hence, to the decay of vortices. 

This work was funded through the New Economy 
Research Fund contract NERF-UOOX0703: Quantum 
Technologies. 
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